【直播】我的基因组82：如何对maf格式的突变文件统计vaf

前面我们已经详细讲解了vaf，而且也对我的vcf文件进行了统计。【直播】我的基因组81：看看我的vcf文件的vaf分布情况

文末我提到了一篇文章里面画出了TCGA的一个tumor样本的somatic mutation的vaf分布！

后台很多小伙伴留言那个该如何做。

其实非常简单的。

前面我们介绍了TCGA的maf文件的下载。

我们随便看一个癌种吧！

先在shell里面

1. cut -f 40-42  TCGA.BLCA.mutect.somatic.maf >TCGA.BLCA.counts

然后在R里面：

1. a=read.table("TCGA.BLCA.counts",header = T,stringsAsFactors = F)

2. a$vaf=a$t_alt_count/a$t_depth

3. hist(a$vaf,breaks = 0:100/100)

就出图啦：

虽然很粗糙，但是意思已经明白了。

但是可以看出tumor里面的vaf分布其实已经不再是扔硬币那样的概率了，对于杂合位点来说。

原因很多，首先tumor不一定是单纯的二倍体了，其次tumor样品一般来说本身异质性高，而我们测序是混合多个细胞的，有一些突变有一些并不突变。

而且纯合的somatic mutation几乎没有，因为somatic mutation是tumor过滤了normal后留下来的变异位点，不是遗传多样性，突变这个过程既然是后天产生的，就很难保证取样部分的几百万个细胞全部突变了。

为什么取对TCGA下载的maf文件取第40~42列呢，也是经过搜索学习的。

MAF格式定义如下：

但是它只规定了前34列，后面所有的列，都是可以自由扩展的。

而我们下载的maf文件共有114 columns，如下：

Hugo_Symbol

Entrez_Gene_Id

Center

NCBI_Build

Chromosome

Start_Position

End_Position

Strand

Variant_Classification

Variant_Type

Reference_Allele

Tumor_Seq_Allele1

Tumor_Seq_Allele2

dbSNP_RS

dbSNP_Val_Status

Tumor_Sample_Barcode

Matched_Norm_Sample_Barcode

Match_Norm_Seq_Allele1

Match_Norm_Seq_Allele2

Tumor_Validation_Allele1

Tumor_Validation_Allele2

Match_Norm_Validation_Allele1

Match_Norm_Validation_Allele2

Verification_Status

Validation_Status

Mutation_Status

Sequencing_Phase

Sequence_Source

Validation_Method

Score

BAM_File

Sequencer

Tumor_Sample_UUID

Matched_Norm_Sample_UUID

HGVSc

HGVSp

HGVSp_Short

Transcript_ID

Exon_Number

t_depth

t_ref_count

t_alt_count

n_depth

n_ref_count

n_alt_count

all_effects

Allele

Gene

Feature

Feature_type

Consequence

cDNA_position

CDS_position

Protein_position

Amino_acids

Codons

Existing_variation

ALLELE_NUM

DISTANCE

TRANSCRIPT_STRAND

SYMBOL

SYMBOL_SOURCE

HGNC_ID

BIOTYPE

CANONICAL

CCDS

ENSP

SWISSPROT

TREMBL

UNIPARC

RefSeq

SIFT

PolyPhen

EXON

INTRON

DOMAINS

GMAF

AFR_MAF

AMR_MAF

ASN_MAF

EAS_MAF

EUR_MAF

SAS_MAF

AA_MAF

EA_MAF

CLIN_SIG

SOMATIC

PUBMED

MOTIF_NAME

MOTIF_POS

HIGH_INF_POS

MOTIF_SCORE_CHANGE

IMPACT

PICK

VARIANT_CLASS

TSL

HGVS_OFFSET

PHENO

MINIMISED

ExAC_AF

ExAC_AF_AFR

ExAC_AF_AMR

ExAC_AF_EAS

ExAC_AF_FIN

ExAC_AF_NFE

ExAC_AF_OTH

ExAC_AF_SAS

GENE_PHENO

FILTER

src_vcf_id

tumor_bam_uuid

normal_bam_uuid

GDC_Validation_Status

GDC_Valid_Somatic

而能够统计vaf的列，我已经标记了红色，就是40~42列咯。